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Using Green-function many-body theory, we present the details of a formally exact adiabatic-connection 
fluctuation-dissipation density-functional theory based on range separation, which was sketched in Toulouse, 
Gerber, Jansen, Savin and Angyan, Phys. Rev. Lett. 102, 096404 (2009). Range-separated density-functional 
theory approaches combining short-range density functional approximations with long-range random phase ap- 
proximations (RPA) are then obtained as well-identified approximations on the long-range Green-function self- 
energy. Range-separated RPA-type schemes with or without long-range Hartree-Fock exchange response kernel 
are assessed on rare-gas and alkaline-earth dimers, and compared to range-separated second-order perturbation 
theory and range-separated coupled-cluster theory. 



I. INTRODUCTION 

Range-separated density-functional theory has emerged as 
a powerful approach for improving the accuracy of standard 
Kohn-Sham (KS) density-functional theory [1,2] applied with 
usual local or semi-local density-functional approximations, 
in particular for electronic systems with strong (static) or 
weak (van der Waals) correlation effects. Based on a sepa- 
ration of the electron-electron interaction into long-range and 
short-range components, it permits a rigorous combination of 
a long-range explicit many-body approximation with a short- 
range density-functional approximation (see, e.g., Ref. 3 and 
references therein). Several many-body approximations have 
been considered for the long-range part: configuration in- 
teraction [4, 5], multi-configuration self-consistent-field the- 
ory [6-8], second-order perturbation theory [9-13], coupled- 
cluster theory [14—18], multi-reference second-order pertur- 
bation theory [19], and several variants of the random phase 
approximation (RPA) [20-24]. 

In the context of the recent revived interest in RPA-type ap- 
proaches to the electron correlation problem in atomic, molec- 
ular and solid-state systems [25-48], several range-separated 
approaches using long-range RPA-type approximations have 
indeed been proposed and show promising results, in particu- 
lar for describing weak intermolecular interactions. Toulouse 
et al. [20] have presented a range-separated RPA-type the- 
ory including the long-range Hartree-Fock exchange response 
kernel. Janesko et al. [21-23] have proposed a simpler range- 
separated RPA scheme with no exchange kernel and in which 
the RPA correlation energy has been rescaled by an empir- 
ical coefficient. Paier et al. [24] have added the so-called 
second-order screened exchange to the latter scheme, which 
appears to correct the self-interaction error. In all these cases, 
range separation tends to improve the corresponding full- 
range RPA-type approach, avoiding the inaccurate description 
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and slow basis-set convergence of short-range correlations in 
RPA. 

In Ref. 20, only the main lines of range-separated density- 
functional theory with long-range RPA were presented. In 
this work, we give now all the missing details of the the- 
ory. Using Green-function many-body theory, we construct 
a formally exact adiabatic-connection fluctuation-dissipation 
density-functional theory based on range separation, with- 
out the need of maintaining the one-particle density constant. 
Range-separated RPA-type schemes are then obtained as well- 
identified approximations on the long-range Green-function 
self-energy. The range-separated RPA-type methods with or 
without long-range Hartree-Fock exchange response kernel 
are assessed on rare-gas and alkaline-earth dimers, and com- 
pared to range-separated second-order perturbation theory and 
range-separated coupled-cluster theory. The most tedious de- 
tails of the theory are given in the appendices. 



II. THEORY 

A. Range-separated density-functional theory 

In range-separated density-functional theory (see, e.g., 
Ref. 3), the exact ground-state energy of an A^-electron system 
is expressed as the following minimization over multidetermi- 
nant wave functions 

£ = imn [{^\f + V„, + Wl\'i') + E^Jn^]] , (1) 

where T is the kinetic energy operator, y„(, is 
the nuclei-electron interaction operator, W^^ - 
(1/2) jJ"(iriiir2W^'"^(ri2)«2(ri,r2) is a long-range (k) 
electron-electron interaction written with w^l^{r) - erf(//r)/r 
and the pair-density operator «2(ri,r2), and ^^^.^[n] is 
the corresponding //-dependent short-range (sr) Hartree- 
exchange-correlation (Hxc) density functional that Eq. (1) 
defines. The parameter yU in the error function controls the 
range of the separation. The minimizing wave function, 
denoted by ^F'^ yields the exact density. Several approxi- 
mations [3, 7, 14, 18, 49-51] have been proposed for the 
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short-range exchange-correlation (xc) functional £'j^.[m], and 
an approximate scheme must be used for the long-range wave 
function part of the calculation. 

In a first step, the minimization in Eq. (1) is restricted to 
single-determinant wave functions <1>, leading to the range- 
separated hybrid (RSH) approximation [9] 

£rsh = mm {<0|f + y„, + W^jO) + E'^Jn^]} , (2) 

which does not include long-range correlation. The mini- 
mizing determinant (Do is given by the self-consistent Euler- 
Lagrange equation 

H^m = £ol<I>o), (3) 

where £o is the Lagrange multiplier for the normalization con- 
straint and Ho is the RSH reference Hamiltonian 

A) = f + V„e + Vh.v,hf[1'o] + Vh,.[«Oo], (4) 

which includes the Hartree-Fock (HF)-type long-range 
Hartree-exchange (Hx) potential V^\.j]p[Oo] and the short- 
range local Hxc potential V'^'^^.p[n] = Jc/rvy^.^[n](r)«(r) 
written with v^^^[n](r) = 5£'^^.^[«]/&i(r) and the den- 
sity operator «(r). As usual, VhcHF '■^^ 
of a local Hartree part Vjj = J^friv[^(ri)n(ri) with 
^H^*"!) ^ / ^r2wJl"e(ri2)<Oo|«(r2)l<l>o>, and a non-local ex- 
change part V''jjp = ^c/xi(ix2v'j':(xi,X2)«i(x2,Xi) written 

with v';"(xi,X2) = -Wg'X'"i2)<<l'o|ni(xi,X2)|<l>o) and the one- 
particle density-matrix operator ni(xi,X2) expressed with 
space-spin coordinates Xi - (ri, S[) and X2 = (r2, 52)- 

The RSH scheme does not yield the exact energy and den- 
sity, even with the exact short-range functional E'^^^[n]. Nev- 
ertheless, the RSH approximation can be used as a reference 
to express the exact energy as 

E = £rsh + (5) 

defining the long-range correlation energy E^^, for which we 
will now give an adiabatic connection formula. We introduce 
the following energy expression with a formal coupling con- 
stant A 

E, = mm{<^|f + V„, + yJl,,Hp[<l>o] + AW'^m 

+El.M'v]}, (6) 

where the minimization is done over multideterminant wave 
functions *P, TV'' is the long-range M0ller-Plesset-type fluctu- 
ation perturbation operator 

W''^Wl-V'^,^^^[^o], (7) 

and E^^^ is the previously-defined /i-independent short-range 
Hxc functional. The minimizing wave function, denoted by 
^P^J , is given by the self-consistent Euler-Lagrange equation 

HaI'V'^} = £!|I^I[>, (8) 



where is the Lagrange multiplier for the normalization con- 
straint and is the long-range interacting effective Hamilto- 
nian along the adiabatic connection 

= i/„ + ^w'' + (y-,K.]-y-j«o„]). (9) 

For /I = 1, Eq. (6) reduces to Eq. (1), and so the physical 
energy E - E,i=i and density are recovered. For A = 0, the 
minimizing wave function is the RSH determinant ^\^o - <to 
and the Hamiltonian of Eq. (9) reduces to the RSH reference 
Hamiltonian, //^"Lq = //()■ Note that, because the density at 
/I = is not exact, the density necessarily varies along this 
adiabatic connection. Taking the derivative of E,i with respect 
to A, noting that E^ is stationary with respect to "^Jj , and rein- 
tegrating between A = and A = I gives 

£-£,=0+ r dA {¥i\w''\¥^}, (10) 

Jo 

with £,=0 = Cfolf + V'«. + VL,hf[<»o]|1>o)+£^^,J«*o] = Ersh- 
(<l>o|W'' |<l>o). Thus, the long-range correlation energy is 

4'- ^ j^^ dA [{¥^\W''\¥^} - <cDo| W''|cDo>] , (11) 
or, equivalently, 

E^^ - — J" dA J" dxidx2dx\dx'2W^'^(xi,X2;x\,x'2) 

xP|,;,(xi,X2;x;,x;), (12) 

where w''"(xi,X2;Xj,x;!,) - w^^],{ri2)S{xi - x[)S(x2 - x'^) - 
l/iN - 1) [v;^(ri)5(x, - x'l) + v>(x,, x'l)] S{X2 - x^) is the po- 
tential corresponding to the perturbation operator W'' and 
Pj; j(xi,X2;Xj,X2) is the correlation part of the two-particle 
density matrix along the adiabatic connection. 

B. Long-range many-body perturbation theory 

We now derive a formally exact many-body perturbation 
theory to calculate the long-range correlation two-particle 
density matrix Pj,'^. Details are given in Appendix A. 

The one-particle Green function Gj[(l,2) along the adia- 
batic connection of Eq. (9) in terms of space-spin-time co- 
ordinates 1 - (xi,fi) and 2 = (X2,f2) satisfies the following 
Dyson equation 

' (1 , 2) = Go ' ( 1 , 2) - E^'d , 2) - A2-(l , 2), (13) 

where G()(l,2) is the reference Green function correspond- 
ing to the RSH Hamiltonian Hq, SJ[(1,2) is the self-energy 
corresponding to the long-range perturbation operator AW'^ 
and ASY(1,2) is the self-energy correction associated with 
the short-range potential variation term yH j,, ] ~ ^h ic["'I'o] 
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due to the variation of the density [52]. The long-range self- 
energy corresponding to the perturbation operator /ICW^^^ - 
^Hi Hf[*^o]) is decomposed into Hartree, exchange and cor- 
relation contributions as 



S'[(l,2) 



4.c.i[G!f](l,2). 



[Gold, 2) 



= /({2l;,[G\^](l,2)-4jGo](l,2)} 

+5:|:,[g!j](i,2), 



(14) 



where Sj^^[G](l,2) is the sum of a long-range Hartree self- 
energy 

Sl^[G](l,2) = -i J t/3t/4w|;,(l, 3)5(1, 2)5(3, 4)G(4,3+) 
= -iS(l,2) J t/3w;;(l,3)G(3,3+) 

- 6(1,2) j dr^wlimMn) 

= 5(l,2)v|;[«](ri), (15) 

with the instantaneous electron-electron interaction 
Wgj,(l,3) = 5(fi - ?3)w|.g(ri3) and the density extracted 
from the Green function n(r3) - -i'Zjs,G{3,3^) (where 3^ 
stands for x^t^ with = h + tj and tj is an infinitesimal 
positive shift), and a long-range exchange self-energy 

2!,'[G](1,2) = i ^ £/3£/4w;;,(l, 3)5(1, 4)5(2, 3)G(4, 3+) 
= /w|.',(l,2)G(l,2+) 



= -5(fi - f2)w;^(ri2)«i(Xi,X2) 

= 5(fi -f2)v';"[ni](xi,X2), 



(16) 



with the one-particle density matrix extracted from the Green 
function «i(xi,X2) - -/G(xiri,X2f|). The short-range self- 
energy correction corresponding to the operator ^'Hj:t ["'J''i ] ~ 
^Hxc["i'o] written as 

AE-(1,2) = SS^JG![](l,2)-SS.ve[Go](l,2), (17) 

where E^\.JG](1,2) = 5(1, 2)v^.Jn](ri) is the local short- 
range Hxc self-energy. 

The long-range four-point polarization propagator 
x\{^,2\ r,2') along the adiabatic connection is given by the 
solution of the following Bethe-Salpeter-type equation which 
can be derived from the Dyson equation (13) by considering 
variations with respect to Gj[ [see Appendix A, Eq. (A13)] 

{X'()\W\2') = {xl.,y\w\2') 
-Af^^{\,2-V,2') 
-/,';:,( 1, 2; 1', 2'), (18) 

where y,;p ,(1, 2; 1', 2') = -iG^l{\,2')G\i2,V) is an 
independent-particle (IP) polarization propagator, and 
i//;;(l,2; l',2') = /^5SJ^JG|[](1,1')/5G;[(2',2) and 
/^';,(1,2;1',2') = /5S;i;,[G;;](1,1')/5G'J(2',2) are long- 
range Hartree-exchange and correlation kernels. Note that 



these kernels only stem from the self-energy term 2^hj.c^[G'[] 
in Eq. (13) that corresponds to the two-electron interaction 
AW^g, the other self-energy contributions which come from 
the one-electron terms are absorbed in the definition of 
x\{l,2; r, 2'). The Hartree kernel is obtained from Eq. (15) 



/^■( 1, 2; 1 ', 2') = w!.',( 1,2)5(1,105(2, 2') 
= wl>i2)5(fi-f2)5(l,l')5(2,2'), 



(19) 



while the HF-like exchange kernel is obtained from Eq. (16) 

/;'(1,2; l',2') = -w^Ul,2)5(l,2')5(l',2) 

= -wtinim, - f2)5(l,2')5(l',2). (20) 

The fluctuation-dissipation theorem is then used to express 
Pj!^ J as [see Appendix A, Eq. (A24)] 

^e-"'[;rl;"(xi,X2;x;,x^;o.) 
-Xq(^\, X2; x'l , X2; w)]+A{[(xi , X2; x'l , Xj), (21) 

where x^a^^\ ' ^2; x'j , Xj; w) is the frequency-dependent Fourier 
transform of the one-time-interval polarization propagator 
;tJ[(xi,X2;x'j,X2;T = fi - h) = ;\fJ[(xifi,X2f2;Xjf|,x;f+), 
;ifo(x 1 , X2 ; X J , Xj ; w) is the equivalent quantity for the RSH ref- 
erence Hamiltonian (at A - 0), and A'[(xi,X2;Xj,X2) is the 
contribution coming from the variation of the one-particle 
density matrix along the adiabatic connection. The expression 
of Ajj in terms of the Green functions Gj[ and Go is straightfor- 
ward but it is sufficient to write it as Aj[ - r[Gjj ] -r[Go] where 
r is a known functional given in Appendix A [Eq. (A22)]. 

So far, the theory is in principle exact. In the following we 
consider two possible approximations. The RPA approxima- 
tion 



(22) 



corresponds to neglecting long-range exchange-correlation in 
all one-electron properties. Indeed, with this approximation, 
one can check that Gjj = Go is a solution of the Dyson equa- 
tion (13), i.e. the Green function remains unchanged along 
the adiabatic connection. It follows that A'[ - 0, /Jj, ^ = and 
;^l^p,(l,2; l',2') = -/Go(l, 2')Go(2, 1') =;ro(l,2; l',2'). Sim- 
ilarly, the RPAx approximation 

^.^=0, (23) 

corresponds to neglecting long-range correlation only in all 
one-electron properties. Again, this approximation implies 
that the Green function remains unchanged along the adiabatic 
connection, i.e. G'[ = Go and it follows that A'[ = 0, fj^ - and 
x\^ A - Xo- As diff'erent terminologies are used in the quan- 
tum chemistry and condensed-matter physics hterature, let us 
stress that what we call RPA here corresponds to a response 
equation (18) with no exchange-correlation kernel (and it is 
also sometimes called linear response time-dependent Hartree 
theory or direct RPA), and what we call RPAx corresponds 
to a response equation with an additional HF-like exchange 
kernel (and it is also sometimes called linear response time- 
dependent Hartree-Fock theory ox full RPA). 
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C. Expressions in an orbital basis 

The RPA or RPAx equations in an orbital basis are derived 
in details in Appendix B. In the basis of RSH spin orbitals, 
the long-range RPA or RPAx correlation energy writes 

= \ dAy^{ib\wl\ai){Vl,Ujb, (24) 

where / and j refer to occupied spin orbitals, and a and b to 
virtual spin orbitals, {ib\w^l^\aj) are the two-electron integrals 
with long-range interaction, and (Pj!^),aj7, are the matrix el- 
ements of the correlation two-particle density matrix. The 
one-electron terms and v'^ in the perturbation operator in 
Eq. (12) do not contribute to Ej!" because of the occupied- 
virtual/occupied-virtual structure of the two-particle density 
matrix in RPA or RPAx. Following the technique proposed by 
Furche [26], Pj,'^ can be obtained as 



P^'.-(A^B:[y^^(M;;r'^(A;[-B'^)"^-l. (25) 



1/2 , 



>l/2 



with M^"" = (A^"" - B![)'^^ (a;[ + B;[) (A^' - B;[)"^ and the or- 
bital rotation Hessians 

+A [{ib\wl\aj) - ^{ib\w'l\ia)\ , (26a) 



= V^b\wl\ii) - ^{ab\wl\ii)\ . (26b) 

where e, are the RSH orbital eigenvalues, and ^ = or ^ = 1 
for RPA and RPAx, respectively. For spin-restricted closed- 
shell calculations, the correlation energy writes in terms of 
spatial orbitals 

E'^ = i f dAy{ib\wi\aj}CF'^,)i,jh, (27) 
2 ^0 



where / and j now refer to occupied spatial orbitals, and a 

- plr 



and b to virtual spatial orbitals, and 'P|y , is the spin-singlet- 



adapted correlation two-particle density matrix obtained as 

l/2/,_iA-l/2/, .,, ,„i,Al/2 



Ipli- _ o 



('A^'B>0"^CM'f''^CA!l-'B:[)' 



(28) 



with iM|j" = (lAjl" - 1b;[)"''(ia;[ + 'b;[)(ia;[ - 'BIj)'^', and 
the singlet orbital rotation Hessians 

('A^--). = (e„ - £i)3ij6„b 

+A [2{ib\w'l\aj) - ^{ib\w'l\ja)\ , (29a) 

(iB||).^_ = A [2{ab\w':,\ij) - ^{ab\wl\ji)\ . (29b) 

Only singlet excitations contribute to Eq. (27), since the two- 
electron integrals involved vanish for triplet excitations. 



In Eq. (25), it is assumed that Aj[ + Bj[ and Aj[ - Bj[ are pos- 
itive definite. In RPA, this is always the case. On the contrary, 
in RPAx, this is not always the case, i.e. instabilities can be 
encountered, and Eq. (25) can fail. In spin-restricted closed- 
shell formalism, one may encounter singlet instabilities in the 
RPAx theory defined here, for example when dissociating a 
bond, but not triplet instabilities since triplet excitations do 
not contribute at all. In practice, singlet instabilities are usu- 
ally not encountered for weakly-interacting closed-shell sys- 
tems. Note that other variants of RPA-type correlation energy 
expressions using a HE exchange response kernel, such as the 
plasmon formula [38, 53, 54] or the equivalent ring coupled- 
cluster-doubles theory [38], require contributions from both 
singlet and triplet excitations, and are thus subject to triplet 
instabilities (e.g. in a system such as Be2). 

Similarly to the notation used in Ref. 20, the range- 
separated method obtained by adding to the RSH energy the 
long-range RPAx correlation energy - 1 in Eqs. (26) 
or Eqs. (29] will be referred to as RSH-i-lrRPAx. For con- 
sistency, the range-separated method obtained by adding 
to the RSH energy the long-range RPA correlation energy 
[^ = in Eqs. (26) or Eqs. (29)] will be referred to as 
RSH+lrRPA, although it is equivalent to the method called 
"LC-wLDA+dRPA" in Refs. 21-24 in the special case of the 
short-range LDA functional. At second order in the electron- 
electron interaction, the RSH+lrRPAx method reduces to the 
range-separated method of Ref. 9 based on long-range second- 
order M0ller-Plesset perturbation theory, to which we will 
refer as RSH+lrMP2. Since RPA approaches can be seen 
as simple approximations to coupled-cluster theory [38], the 
RSH+lrRPA and RSH+lrRPAx methods bear some resem- 
blance to the range-separated method of Ref. 14 where the 
long-range correlation energy is evaluated by coupled-cluster 
theory (with single, double and perturbative triple excita- 
tions), to which we will refer as RSH+lrCCSD(T). 

We note that one can develop long-rang many-body pertur- 
bation theories starting from other references than the RSH 
reference. For example, starting from the usual (approximate) 
Kohn-Sham reference could be appropriate for solid-state sys- 
tems. For the finite systems considered here, RSH is a good 
reference, as confirmed by other authors [23]. 



m. COMPUTATIONAL DETAILS 

All calculations have been performed with a development 
version of MOLPRO 2008 [55], implementing equations (27)- 
(29). We first perform a self-consistent RSH calculation with 
the short-range PBE xc functional of Ref. 14 (this RSH cal- 
culation could also be referred to as "IrHF+srPBE", a nota- 
tion closer to the one used by other authors [14]) and add 
the long-range MP2, RPA, RPAx or CCSD(T) correlation en- 
ergy calculated with RSH orbitals. For RPA or RPAx, the 
/{-integration in Eq. (27) is done by a 7-point Gauss-Legendre 
quadrature [26]. The range separation parameter is taken at 
// = 0.5 bohr"', in agreement with previous studies [56], with- 
out trying to adjust it for each system. To show the depen- 
dence on the orbitals, the full-range RPA calculations have 
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been done with PBE [57] and HF orbitals, which will be de- 
noted by PBE+RPA and HF+RPA, respectively [58]. The 
full-range MP2, RPAx and CCSD(T) calculations have been 
done with HF orbitals, and thus, for notation consistency, will 
be denoted by HF+MP2, HF+RPAx and HF+CCSD(T), re- 
spectively. We use large Dunning basis sets [59-65]. Core 
electrons are kept frozen in all the full-range and range- 
separated MP2, RPA, RPAx and CCSD(T) calculations (i.e. 
only excitations of valence electrons are considered). The ba- 
sis set superposition error (ESSE) is removed by the coun- 
terpoise method. For the alkaline-earth dimers, it has been 
checked than adding diffuse basis functions or core excita- 
tions do not change significantly the results. Extrapolations 
to the complete basis set (CBS) limit have also been consid- 
ered for some systems. For the full-range methods, the stan- 
dard three-point exponential formula for the HF (or KS) ref- 
erence Eiip{n) - EhfCCBS) +Ae^^" with the cardinal number 
n = 3,4,5, and two-point formula for the correlation energy 
Ecin) = £'f.(CBS) + C/«'' with n -4,5 have been used. For the 
range- separated methods, we have also used these two formu- 
las for the RSH reference and the long-range correlation en- 
ergy, even though in this case the dependence on the cardinal 
number would deserve a detailed study. 

For each dimer interaction energy curve, we choose 16 to 
20 intermolecular distances, with denser sampling around the 
equilibrium distance. A third-order polynomial is used for 
interpolation. The hard core radius is taken as the distance 
where the interaction energy is 0, and the equilibrium dis- 
tance and binding energy are from the minimum of the in- 
terpolated interaction energy curve. The harmonic vibrational 
frequency is obtained from the second-order derivative of the 
energy curve at the equilibrium distance. For Ce dispersion 
coefficients, the interaction energy Ei^t is calculated at seven 
extra distances from 30 to 60 bohr, and the coefficient is 
estimated by averaging with the following formula 



Cfi - exp 



i^(ln|£,„,(i?,)l + 61n(i?,)) 



(30) 



V /=1 

similarly to what has been done in Ref. 22. 

IV. APPLICATIONS 
A. Basis set dependence 

The convergence of the equilibrium binding energy of 
Ar2 with respect to the basis set size up to the CBS 
Umit for the full-range methods HF+MP2, PBE+RPA, 
HF+RPA, HF+CCSD(T) and for the range-separated methods 
RSH+lrMP2, RSH+lrRPA, RSH+lrRPAx, RSH+b-CCSD(T) 
is represented in Fig. 1 . Full-range RPA with PBE orbitals has 
a very strong dependence on the basis size, as already noted 
(e.g. Refs. 20, 26). Full-range RPA with HF orbitals has a 
bit weaker basis dependence, similar to full-range HF+MP2, 
HF+RPAx and HF+CCSD(T). All the range-separated meth- 
ods have essentially identical, very favorable basis set conver- 
gence. Since the slow convergence of full-range methods is 
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FIG. 1: (Color online) Basis set dependence of the equilibrium bind- 
ing energy of Ar2 for different full-range and range-separated meth- 
ods, presented as the percentage of the binding energy recovered with 
respect to the CBS limit (aVTZ, aVQZ and aV5Z stand for aug-cc- 
pVTZ, aug-cc-pVQZ and aug-cc-pV5Z, respectively). 



related to the explicit description of short-range correlation, 
it is not surprising that range-separated methods have a faster 
convergence because they leave the description of short-range 
correlation to the short-range density functional. These results 
are consistent with other studies, e.g. Refs. 22, 24. Note that, 
with the aug-cc-pV5Z basis set, all the range-separated meth- 
ods are essentially converged (98-99% of the CBS binding 
energy), therefore we will not use CBS extrapolations in the 
following. However, one should keep in mind that with this 
basis set the full-range methods are not yet fully converged, 
with about 90% of the CBS binding energy. 



B. Rare-gas dimers 

In Fig. 2, the interaction energy curves of Hea, Nea, Ar2 and 
Kr2, obtained with the full-range and range-separated meth- 
ods are compared. As already known, full-range HF+MP2 
underestimates the interaction energy for the smallest sys- 
tems He2 and Ne2, and overestimates it for the largest sys- 
tems Ar2 and Kr2. Full-range PBE+RPA gives an almost dis- 
sociative curve for He2, and largely underestimates the in- 
teraction energy for Ne2, Ar2 and Kr2. Using HF orbitals 
in full-range RPA drastically improves the interaction energy 
curve for He2, and to a least extend for Ne2, but gives less 
binding for Ar2 and Kr2. Full-range HF+RPAx significantly 
improves over full-range HF+RPA, but still gives underesti- 
mated interaction energies. It can be noted that full-range 
HF+RPAx yields interaction energy curves almost identical to 
the full-range HF+MP2 curves for He2 and Ne2, and almost 
identical to the full-range PBE+RPA curves for Ar2 and Kr2. 
Full-range HF+CCSD(T) gives systematically quite accurate 
interaction energies. Quite similarly to full-range HF+MP2, 
the range-separated RSH+lrMP2 underestimates the interac- 
tion energy for He2 and Ne2, and overestimates it for Ar2 
and Kr2. RSH+lrRPA tends to improve over both full-range 
PBE+RPA and HF+RPA but still leads to significantly under- 
estimated interaction energies. RSH+lrRPAx improves over 



6 



S 



a 

Qi 

§ 

i 



0.01 


-0.01 
-0.02 



-0.03 



-0.04 



1 1 


-t ' 
\He2 


■■■■■■■4-.„ 








^ + +-—^^^ 
























Accurate - 






HF+MP2 






PBE+RPA 


+ 




HF+RPA 






HF+RPAx 


-o 




HF+CCSD(T) 


— B 



0.05 



5 6 7 8 9 

Internuclear distance (bohr) 



Accurate 
HF+MP2 

PBE+RPA 
HF+RPA 

HF+RPAx 
HF+CCSD(T) 



6 7 8 9 

Internuclear distance (bolir) 



Accurate 
HF+MP2 

PBE+RPA 
HF+RPA 

HF+RPAx 
Hr+CCSD(T) 



10 




10 




7 8 9 10 11 12 
Internuclear distance (bohr) 



14 



0.1 









0.1 














0.2 




/ M 








V 




0.3 


11 ® «■ 










0.4 


-11 *«■ i 


fi. 




0.5 


W -//^ 




Accurate 




HF+MP2 


0.6 






PBE+RPA + 






HF+RPA 


0.7 




HF+RPAx O 


0.8 






IIF+CCSD(T) — B — 







8 10 12 14 

Internuclear distance (bohr) 



16 



0.01 



-0.01 



-0.02 



-0.03 



-0.04 

















\ ra ® "/^ Accurate - 




\\f RSH+lrMP2 




\ RSH+lrRPA 


--X--- 


V/ RSH+lrRPAx 


o 


RSH+lrCCSD(T) 


— B 



5 6 7 8 9 

Internuclear distance (bohr) 



Accurate 

RSH+lrMP2 
RSH+lrRPA 
RSH+lrRPAx O 
RSH+lrCCSD(T) 



6 7 8 9 

Internuclear distance (bohr) 



Accurate 
RSH+lrMP2 
RSH+lrRPA 
RSH+lrRPAx 
RSH+lrCCSD(T) 



7 8 9 10 11 12 
Internuclear distance (bohr) 



Accurate 

RSH+lrMP2 
RSH+lrRPA 
RSH+lrRPAx O 
RSH+lrCCSD(T) — B- 



8 10 12 14 

Internuclear distance (bohr) 



10 




10 




14 




16 



FIG. 2: (Color online) Interaction energy curves of He2, Ne2, Ar2 and Kr2 calculated by different full-range (left) and range- separated (right) 
methods. The basis is aug-cc-pV5Z. The accurate curves are from Ref. 66. 



7 



TABLE I: Hard-core radii cr (bohr), equilibrium distances (bohr), equilibrium binding energies (mhartree), harmonic vibrational fre- 
quencies a>g (cm"') and dispersion coefficients C(, for ten homonuclear and heteronuclear rare-gas dimers from different full-range and range- 
separated methods with aug-cc-pV5Z basis. Mean absolute percentage errors (MA%E) are also given. 







HF+MP2 


PBE+RPA 


HF+RPA 


HF+RPAx 


HF+CCSD(T) 


RSH+lrMP2 


RSH+lrRPA 


RSH+lrRPAx 


RSH+lrCCSD(T) Estimated exact" 


He2 






















cr 


5.20 


6.81 


5.34 


5.18 


5.03 


5.35 


5.39 


5.25 


5.17 


5.02 


Re 


5.83 


8.16 


5.95 


5.82 


5.65 


6.00 


6.10 


5.92 


5.85 


5.62 


De 


0.0208 


0.0021 


0.0145 


0.0218 


0.0313 


0.0202 


0.0183 


0.0255 


0.0309 


0.0348 


(lie 


26.9 


4.5 


24.1 


27.4 


33.6 


26.2 


22.3 


28.6 


30.4 


34.3 


Cf, 


1.13 


1.36 


0.88 


1.14 


1.46 


1.42 


1.34 


1.67 


1.91 


1.461 


He-Ne 






















tr 


5.32 


5.81 


5.44 


5.29 


5.13 


5.33 


5.38 


5.27 


5.19 


5.16 


«. 


5.95 


6.37 


6.08 


5.91 


5.77 


5.99 


6.07 


5.93 


5.87 


5.76 


£>. 


0.0401 


0.0064 


0.0284 


0.0410 


0.0609 


0.0458 


0.0401 


0.0533 


0.0638 


0.0660 


big 


28.8 


13.0 


23.8 


29.5 


34.3 


28.4 


26.2 


30.9 


33.5 


36.1 


Cf, 


2.43 


2.77 


1.84 


2.32 


3.07 


3.12 


2.84 


3.44 


4.04 


3.029 


He-Ar 






















(T 


6.02 


6.31 


6.27 


6.11 


5.92 


6.01 


6.14 


5.99 


5.87 


5.92 


Re 


6.73 


6.96 


6.97 


6.83 


6.64 


6.77 


6.89 


6.73 


6.63 


6.61 


De 


0.0736 


0.0307 


0.0424 


0.0608 


0.0874 


0.0808 


0.0616 


0.0854 


0.1071 


0.0937 


aj^ 


32.3 


24.1 


25.9 


29.4 


35.7 


31.5 


29.0 


33.3 


37.4 


36.0 


Ce 


9.1 


9.1 


6.1 


7.6 


11.6 


10.6 


8.7 


10.8 


12.6 


9.538 


He-Ki- 






















CT 


6.38 


6.67 


6.67 


6.50 


6.28 


6.35 


6.52 


6.34 


6.22 


6.25 


Re 


7.15 


7.37 


7.42 


7.26 


7.05 


7.14 


7.31 


7.13 


7.03 


6.98 


De 


0.0747 


0.0337 


0.0423 


0.0606 


0.0881 


0.0833 


0.0613 


0.0857 


0.1084 


0.0996 


(lie 


30.1 


22.3 


23.4 


26.3 


32.4 


30.7 


25.9 


31.2 


34.2 


33.7 


C6 


12.9 


12.5 


8.5 


10.7 


14.0 


14.9 


12.0 


14.7 


17.3 


13.40 


Ne2 






















tr 


5.47 


5.63 


5.57 


5.43 


5.28 


5.36 


5.43 


5.33 


5.27 


5.23 


Re 


6.11 


6.18 


6.19 


6.07 


5.90 


6.03 


6.10 


5.98 


5.93 


5.84 


De 


0.079 


0.037 


0.056 


0.077 


0.118 


0.102 


0.088 


0.111 


0.131 


0.134 


(lie 


22.8 


18.7 


19.7 


22.6 


28.8 


23.8 


22.9 


25.9 


28.3 


29.4 


Cf, 


5.24 


6.84 


3.91 


4.77 


6.35 


6.80 


6.10 


7.03 


8.08 


6.383 


Ne-Ai- 






















cr 


6.02 


6.21 


6.28 


6.13 


5.94 


5.92 


6.06 


5.93 


5.84 


5.89 


Re 


6.72 


6.87 


7.01 


6.85 


6.65 


6.66 


6.80 


6.67 


6.59 


6.57 


De 


0.163 


0.095 


0.092 


0.126 


0.189 


0.196 


0.147 


0.192 


0.235 


0.211 


(iJe 


25.3 


21.6 


17.4 


22.6 


27.7 


27.2 


23.0 


26.9 


29.3 


28.7 


C6 


19.2 


18.9 


12.5 


15.2 


18.2 


22.6 


18.3 


21.8 


25.3 


19.50 


Ne-Kr 






















O" 


6.31 


6.53 


6.61 


6.46 


6.24 


6.20 


6.36 


6.23 


6.14 


6.17 


Re 


7.08 


7.21 


7.36 


7.20 


6.98 


6.97 


7.13 


7.01 


6.91 


6.89 


De 


0.174 


0.104 


0.096 


0.131 


0.201 


0.212 


0.153 


0.201 


0.248 


0.224 


Itln 


22.4 


19.0 


17.0 


19.8 


24.5 


24.4 


20.7 


23.1 


26.5 


25.3 


Ci 


27.0 


26.2 


17.4 


21.1 


27.4 


31.5 


24.8 


29.5 


34.0 


27.30 


Al2 






















tr 


6.32 


6.61 


6.74 


6.60 


6.41 


6.32 


6.55 


6.40 


6.28 


6.37 


Re 


7.10 


7.36 


7.52 


7.37 


7.17 


7.11 


7.34 


7.18 


7.07 


7.10 


De 


0.483 


0.269 


0.215 


0.289 


0.414 


0.484 


0.308 


0.420 


0.542 


0.454 


OJe 


32.7 


25.5 


21.4 


25.5 


30.7 


32.1 


25.5 


30.0 


33.5 


32.1 


Cs 


76.3 


58.6 


42.9 


52.0 


64.5 


80.7 


57.4 


69.6 


85.0 


64.30 



Ar-Kr 



tr 
Re 

De 

OJg 
C6 


6.55 
7.36 

0.570 
29.5 

109.9 


6.85 
7.64 
0.319 
22.9 
82.1 


7.00 
7.81 
0.248 
19.4 
60.7 


6.85 
7.66 
0.334 
22.7 
73.6 


6.65 
7.45 
0.481 
27.3 
94.8 


6.55 
7.37 

0.563 
28.7 

114.1 


6.80 
7.62 
0.346 
22.5 
80.0 


6.64 
7.46 
0.472 
26.3 
97.4 


6.52 
7.34 

0.615 
29.8 

117.1 


6.59 
7.35 

0.531 
28.6 

91.13 


Kr2 

(T 
Re 
De 

Ce 


6.77 
7.60 
0.691 
25.1 
159 


7.09 
7.90 
0.388 
19.8 
116 


7.24 
8.08 
0.296 
16.2 
86 


7.10 

7.92 
0.396 
19.7 
105 


6.88 
7.70 
0.575 
23.2 
132 


6.77 
7.61 
0.671 
24.4 
162 


7.05 
7.89 
0.397 
19.2 
109 


6.88 
7.72 
0.542 
21.9 
134 


6.75 
7.60 
0.713 
25.0 
163 


6.79 
7.58 

0.638 
24.4 

129.6 


MA%E (%) 

tr 

Re 

De 

Ue 

Ce 


2.1 
2.1 
23 
12 
13 


9.3 
9.4 
62 
36 
7.0 


6.3 
6.1 
56 
33 
36 


3.8 
3.9 
39 
21 
22 


0.7 
1.0 
10 
3.4 
4.1 


1.8 
2.2 

16 
9.5 

14 


3.9 
4.5 
36 
23 
9.2 


1.5 
2.2 
14 
10 
10 


1.0 
1.0 
11 
4.6 

29 


0.0 
0.0 
0.0 
0.0 
0.0 





" From Ref . 66 



both RSH-i-lrRPA and full-range HF+RPAx; it still systemati- 
cally underestimates the interaction energy at equilibrium, but 
appears quite accurate at medium and large distances. On the 
contrary, RSH-i-lrCCSD(T) systematically overestimates the 



interaction energy at medium and large distances. 

The hard-core radii, equilibrium distances, equilibrium 
binding energies, harmonic vibrational frequencies and dis- 
persion coefficients Ce for ten homonuclear and heteronu- 
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FIG. 3; (Color online) Interaction energy curves of Be,, Mg, and Ca2 calculated by full-range (left) and range-separated (right) methods. The 
basis is cc-pV5Z. The accurate curves are from Refs. 67, 68 and 69. 



clear rare-gas dimers calculated with the full-range and range- 
separated methods are given in Table I. The trends seen in 
Fig. 2 are confirmed. Full-range RPA (with PBE or HF or- 
bitals) yields very inaccurate equilibrium properties. Full- 
range HF+RPAx improves over full-range HF-)-RPA (with 
the exception of C(, coefficients which turn out to be quite 
good in PBE+RPA for these systems) but the errors remain 
large. Range separation largely improves RPA and RPAx. 
RSH+lrRPAx gives much better equilibrium properties than 
RSH+lrRPA, with mean absolute percentage errors smaller 
by more than a factor of two, while these two methods give 
similar accuracy for Ce coefficients. Full-range HF+MP2 is 
reasonably accurate and range separation has a much smaller 
impact on it. For these systems, RSH+lrMP2 gives an over- 



all similar accuracy than RSH+RPAx, although the Cg coef- 
ficients tend to be globally more accurate in RSH+lrRPAx. 
Full-range HF+CCSD(T) gives the best results. Surpris- 
ingly, range separation tends to deteriorate the accuracy of 
CCSD(T), especially forCg coefficients. Nevertheless, among 
the range-separated methods, RSH+lrCCSD(T) still gives the 
best equilibrium properties. 



C. Alkaline-earth dimers 

In Fig. 3, the interaction energy curves of Be2, Mg2 
and Ca2, obtained with the full-range and range-separated 
methods are compared. These systems have static corre- 
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TABLE II: Hard-core radii cr (bolir), equilibrium distances (bohr), equilibrium binding energies £)<, (mhartree), harmonic vibrational 
frequencies cj^ (cm"' ) and dispersion coefficients Ct for Be2, Mgi and Ca.2 from different full-range and range-separated methods with cc-pV5Z 
basis. Mean absolute percentage errors (MA%E) are also given. 
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lation effects, especially Be2, and are thus more challeng- 
ing for the single-reference methods tested here. Full-range 
PBE-hRPA gives unphysical interaction energy curves, with 
a large bump for Be2, and with essentially no bond for 
Mg2 and Ca2. Full-range HF-i-RPA yields an almost dis- 
sociative curve for Be2 with no bump (which is consistent 
with Ref. 43), and physically reasonable curves for Mg2 and 
Ca2. Full -range HF-hRPAx moderately improves over full- 
range HF+RPA. Among the full-range methods, HF+MP2 
and HF+CCSD(T) clearly give the best interaction energy 
curves. As for rare-gas dimers, RSH+lrRPA always largely 
underestimates the interaction energy. RSH+lrMP2 and 
RSH+lrRPAx give much less underestimated interaction en- 
ergies, with RSH-i-lrMP2 being a bit more accurate for Mg2 
and Ca2. While RSH+lrCCSD(T) largely overestimates the 
interaction energy for Be2, it is remarkably accurate for Mg2 
and Ca2. We note that RSH+lrCCSD(T) could be made more 
accurate for Be2 by choosing a larger range-separation param- 
eter [71]. 



The hard-core radii, equilibrium distances, equilibrium 
binding energies, harmonic vibrational frequencies and dis- 
persion coefficients Ce for Be2, Mg2 and Ca2 are given in 
Table II. It is confirmed that range separation largely im- 
proves the equilibrium properties of RPA and RPAx. Again, 
RSH+lrRPAx is much more accurate than RSH+lrRPA, with 
mean absolute percentage errors smaller by about a factor of 
two. Range separation also overall brings a significant im- 
provement in MP2. Among the range-separated methods, 
RSH+lrCCSD(T) gives the best equilibrium properties. 



V. CONCLUSIONS 

We have expounded the details of a formally exact 
adiabatic-connection fluctuation-dissipation density- 
functional theory based on range separation. Range-separated 
density-functional theory with random phase approximations 
including or not the long-range Hartree-Fock exchange re- 
sponse kernel (referred to as RSH+lrRPA and RSH+lrRPAx, 
respectively) are then obtained as well-identified approx- 
imations on the long-range Green-function self-energy 
[Eqs. (22) and (23)]. The long-range Green function does not 
vary along the adiabatic connection at the RSH+lrRPA and 
RSH+lrRPAx levels, which makes these schemes relatively 
simple compared to the exact theory. In practice, RSH+lrRPA 
and RSH+lrRPAx have been applied in a spin-restricted 
closed-shell formalism, in which both schemes only include 
spin-singlet orbital excitations, and thus are not subject to 
triplet instabilities. 

These range-separated RPA-type schemes have been tested 
on rare-gas and alkaline-earth dimers, featuring challenging 
weak (van der Waals) interactions. Both range separation and 
inclusion of the exact Hartree-Fock response kernel largely 
improve the accuracy of RPA. The RSH+lrRPAx method ap- 
pears as a reasonably accurate method for weak interactions, 
but globally less accurate for equilibrium properties than 
the more intensive range-separated coupled-cluster method. 
Although, for the small systems considered here, range- 
separated second-order perturbation theory (RSH+IrMP2) 
turns out to yield results similarly as accurate as those from 
RSH+lrRPAx (and in fact more accurate for Mg2 and Ca2), 
a recent investigation [72] shows that RSH+lrRPAx corrects 
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the overestimation of the binding energy in RSH+lrMP2 for 
larger weakly-interacting stacked complexes, such as the ben- 
zene dimer 
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Appendix A: Adiabatic-connection fluctuation-dissipation 
density-functional theory 

In this appendix, we outline a general, formally 
exact adiabatic-connection fluctuation-dissipation density- 
functional theory, using Green-function many-body theory. 
For further details on standard Green's function theory, see 
e.g. Refs. 73-76. 



1. Adiabatic connection 

We consider the following adiabatic connection defined by 
the /i-dependent energy 

Ea = mm [{'V\k(i + AW\^') + F[n^,]] , (Al) 

where K() is an arbitrary one-particle Hamiltonian, W is a 
perturbation operator (generally, the sum of a two-particle 
operator Wee and an one-particle operator) and F[n] is a A- 
independent density functional. The minimizing multideter- 
minant wave function satisfies the Euler-Lagrange equa- 
tion 

HA\'i'A)^&AWA), (A2) 

where £^ is the Lagrange multiplier for the normalization con- 
straint, and Hx is the effective Hamiltonian along the adiabatic 
connection 

Ha^Kq + AW + Va, (A3) 

where Va - J dr6F[n'^J/dn(r)h(r) is a self-consistent one- 
particle potential operator Note that Ha=i is not necessar- 
ily the physical Hamiltonian. This adiabatic connection links 
the energy of interest £^=1 to the reference energy £^=0 - 
(OqIA'oIOo) + F[nij,J calculated with the single-determinant 
wave function Oq = *l'^=o of the reference Hamiltonian = 
Kq + V(). The one-particle density is not kept constant with 
respect to A. 

An adiabatic connection formula for Ea=] is found by tak- 
ing the derivative of Ea with respect to A, noting that Ea is 



stationary with respect to and reintegrating between A - 
and A = 1 

EA=i=EA=o+ f dA{^'A\W\^>A). (A4) 
Jo 

The correlation energy, defined as Ec - Ea=\ - Ea=o - 
(dEA/dA)A=o where {dEA/dA)A=() = {<i>o\W\^o) is the first- 
order energy correction, is thus given by 

Ec^ dA [(T.IWI*!',,) - (Oo|W|cDo)] , (A5) 

or, equivalently, in the representation of space-spin coordi- 
nates X - (r, s) 

Ec - — J" dA J" dxidx2dx\dx'2w{xi,X2;x\,x'2) 

XPe,A{Xi,X2;x\,x'2), (A6) 

where w(xi,X2;Xj,X2) is the interaction potential correspond- 
ing to the operator W and f X2; x'j, x^,) is the correlation 
part of the two-particle density matrix along the adiabatic con- 
nection. 

This exposition encompasses both standard full-range 
many-body theory and range-separated density-functional 
theory. Indeed, if Kq is the Hartree-Fock Hamiltonian (i.e., 
ko - f + V„e + Vha,hf), W is the standard M0ller-Plesset 
fluctuation perturbation operator (i.e., W - Wee - V'ha.hf) and 
F[n] = then Eq. (A6) yields the full-range many-body cor- 
relation energy, defined with respect to the Hartree-Fock en- 
ergy. Similarly, with the corresponding long-range operators 
ko = t + V„e + y^^^ HF and W = Wi - y^'^. HP and the short- 
range density functional F[n] - E^^^[n], Eq. (A6) yields now 
the long-range correlation energy, defined with respect to the 
RSH energy [Eq. (5)]. 

2. One-particle Green function 

The one-particle Green function along the adiabatic con- 
nection is defined as 

G,,(i,2) = -a^Ainhm^Ai^JWA), (A7) 

where 1 - (xi,fi) and 2 - (X2,f2) refer to space-spin and 
time coordinates, 0-,i(l) - e'"-'''ifi{xi)e''"^'' and ip\(2) - 
g'H,i'2^t(-x2)e"'^''- are the annihilation and creation operators 
in the Heisenberg picture, and T is the Wick time-ordering 
operator 

A Dyson-type equation connects the inverse of Ga to the 
inverse of the Green function associated with the one-electron 
Hamiltonian Kq + Va, denoted by Gv,a, 

G:,\l,2)^Gy\(l,2)-^A(h2), (A8) 

which can be considered as the definition of the self-energy 
Z ). In turn, the inverse of Gv,a can be expressed from the in- 
verse of the Green function Go of the reference Hamiltonian 
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Ho - Ko + V() as Gy\ = G,,' - [v,i - vq], where and vq are 
the one-electron potentials associated with y,i and Vq, respec- 
tively. 

For time-independent Hamiltonians, the Green function 
only depends on the time difference t = ti - t2, so one de- 
fines G^(xi,X2;t) = G^(xifi,X2f2), which has a discontinu- 
ity at T = 0. The one-particle density matrix «i_^(xi,X2) = 
(*l'^|«i(xi,X2)|*l'^), with «i(xi,X2) = i^"''(x2)iA(xi), can be ob- 
tained from the limit r — » 0" 

ni,i(xi,X2) = -/G^(xi,X2;t = 0"). (A9) 



3. Four-point polarization propagator 

The four-point polarization propagator along the adiabatic 
connection is defined as 

Xa(1,2; 1',2') = /[G2,,(l,2; l',2') - G,,(l, l')G,(2,2')] , 

(AlO) 

where G2,,i is the two-particle Green function 

G2,,(l,2; l',2') = ~{^AT[hWh(m\(2')'lr\{l')WA}, 

(All) 

Alternatively, using the Schwinger derivative technique, xa 
can be expressed as the functional derivative of the one- 
particle Green function with respect to the two-point potential 
x'A (see, e.g., Refs. 73, 76) 



;^,(1,2;1',2') = -/ 



^Gv.,i(l,r) 



(A12) 



The four-point polarization propagator satisfies a so-called 
Bethe-Salpeter equation that directly stems from the Dyson 
equation of Eq. (A8). Considering variations with respect to 
iGa (achieved through variations of v^) yields 



6G-\1,1') SGy\(l,l') .51,(1,1') 



+ ;- 



5Gi(2',2) 5Gi(2',2) 5Ga{2',2) 



(A13) 



The term on the left-hand side of Eq. (A 13) gives straightfor- 
wardly 



^^7(1,1') , , , , 

.J = iG:^'(i,2')G:^'(2,i') 



SGa{2',2) 



= ;^-p',,(l,2;l',2'), 



(A14) 



where ;^'iP,i(l, 2; 1', 2') = -iG^(l, 2')G^(2, 1') is a so-called 
independent-particle (IP) polarization propagator [77]. The 
first term on the right-hand side of Eq. (A 13) gives the in- 
verse of the four-point polarization propagator, according to 
Eq. (A12), 

-^liiM - '■^G:K2^2^-^"(^'2^^'2>' (^1^^ 

and the second term is the so-called Bethe-Salpeter four-point 
kernel 



52,(1, 1') , , 



(A16) 



and finally, using Eqs. (A14)-(A16) in Eq. (A13), the Bethe- 
Salpeter equation for;\^, writes 

;rl'(l,2;l',2')=A'rp„i(l,2;l',2')-/^(l,2;l',2'). (A17) 



4. Fluctuation-dissipation theorem 

Similarly to the expression of the one-particle density ma- 
trix in terms of the one-particle Green function [Eq. (A9)], 
the two-particle density matrix can be extracted from the 
polarization propagator. Defining ;i',(xi,X2;Xj,X2;T) = 
;ifi(xifi , X2r2; X;?!, XjfJ), i.e. the polarization propagator with 
times fj — > f| and fj f J which depends only on the time dif- 
ference T - t[ - t2, it is easy to check that in the limit t — > 0", 
after applying the time-ordering operator in Eq. (All) and us- 
ing Eq. (A9), one has the following relation 

;>f,i(xi,X2;x'i,X2;r = 0") = <*P,)|«i(x2,X2)ni(xi,x'i)|'P^) 

-ni,^(xi,x'i)ni,,(x2,X2). 

(A18) 

The two-particle density matrix n2,.i(xi,X2;x'j,X2) = 
(^^|0'"'"(x2)l^"'"(x^)0'(xl)0•(x2)|*P.^) can thus be expressed as 

n2„i(xi,X2;x'i,X2) = (*l'^|ni(x2,X2)«i(xi,x'i)|1',) 
-6{x\ -X2)ni,^(xi,X2) 
= i>^(xi,X2;x'i,X2;T = 0") 

+ni,^(Xi,x'i)«i,,(X2,X2) 

-5(x;-X2)«i,,(xi,x^). (A19) 

The correlation part of the two-particle density matrix P^j = 
n2,A - n2,A=o is thus 

f'c„t(xi,X2;Xi,X2) = !>',)(xi,X2;x'i,X2;t = 0") 
-i>fo(xi,X2;x'i,X2;T = 0") 
+A,(xi,X2;x'i,x^), (A20) 

where xo is the polarization propagator of the non-interacting 
reference system for A - Q, and is a term coming from the 
variation of the one-particle density matrix along the adiabatic 
connection 

A^(xi,X2;Xj,X2) = ni,^(xi,Xj)ni,^(x2,X2) 
-S{x\ -X2)ni.,i(xi,X2) 
-ni,o(xi,x'i)ni,o(x2,X2) 
+5(x;-X2)ni,o(xi,x^). (A21) 

Using Eq. (A9), one can also express this term with the Green 
function as A ) = r[G,)]-r[Go] where we define the functional 

r as 

r[G] = -G(xi,x;;t = 01G(x2,X2;t = 0") 

+d{x\ - X2)!G(xi, x^; T = 0"). (A22) 

Finally, introducing the Fourier transform of 
;if^(xi , X2; X J , Xj; r) in terms of the frequency cj, 

oc 2m 

x;ifi(xi,X2;Xi,X2; w), (A23) 
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we arrive at the form of the fluctuation-dissipation that we use 

2^. e'"^* (x 1 , X2 ; x'l , X2 ; cu) 
o(x 1 , X2 ; X 1 , X2 ; w)] + (x 1 , X2 ; X 1 , Xj) . ( A24) 

Appendix B: Random phase approximation in an orbital basis 

In this appendix, we give the working equations in an or- 
bital basis resulting from the many-body theory outlined in 
Appendix A, in the special case of a random phase approxi- 
mation (RPA)-type simplification. For further details, see e.g. 
Refs. 26, 53, 78, 79. 

1. Expressions in a spin-orbital basis 

In the RPA and RPAx approximations, the Green function 
does not vary along the adiabatic connection, i.e. - Go, 
which imphes that the independent-particle polarization prop- 
agator [Eq. (A14)] is just the non-interacting reference polar- 
ization propagator, ;ifiP,^(l, 2; 1', 2') = -;Go(l, 2')Go(2, 1') - 
Xo{l,2; l',2'), and in the fluctuation-dissipation theorem of 
Eq. (A24) the term coming from the variation of the one- 
particle density matrix vanishes, A^ = 0. 

The frequency-dependent non-interacting polarization 
propagator has the following well-known Lehmann represen- 
tation 



(xi,X2;x'i,X2; w) = ^ 



a;(x;)0„(xi)0:(x^)^,(x2) 



0*(x^)0„(x2)C(x'i)0-(xi) 
w + (e,, - 6,) - ;0+ 



(Bl) 



where (ppix) and Cp are the spin orbitals and corresponding 
eigenvalues of the reference system, and / and a run over oc- 
cupied and virtual spin orbitals, respectively. Hence, xo can 
be completely represented in the basis of spin-orbital prod- 
ucts, (f)*,{x[)<pg{xi), where p refer to an occupied orbital and q 
to a virtual orbital, and vice versa, with matrix elements 

(ino(w))pf^,„ = J dxidx2dx[dx'2(f)p{x\)(p*g(xi) 

XXoi^l . X2; X'l , X2; W)0;(X2)0,(X2). (B2) 

Assuming orthonormality of the spin orbitals, the matrix ele- 
ments are easily calculated 



Ll> - (€„ - e,) + /0+ 



5n5. 



ijOab 



a> + (e„ - e,) - /0+ 



(ino(w))a,jfc = (ino(w)),v,_ij = 0, 



(B3a) 



(B3b) 



(B3c) 



where both / and j refer to occupied orbitals and both a and b 
to virtual orbitals. The matrix is thus diagonal, and the inverse 
of xo has the following 2x2 supermatrix representation 



Ae 
Ae 



1 

-1 



(B4) 



where Aeiajb - (fo - edSijdab, each block matrices being re- 
indexed with the composite indices ia and jb. 

In the RPA and RPAx approximations, the Bethe-Salpeter 
kernel of Eq. (A16) is approximated as the frequency- 
independent Hartree(-Fock) form [Eqs. (19) and (20)] 

fA{XuX2\x[,x'2) = AWee(ri2)[6{Xi - X[)6(X2 - X'^) 

-^5(xi-x^)5(x;-X2)], (B5) 

where Wedrn) is a two-particle interaction, and ^ = or ^ = 1 
for RPA and RPAx, respectively. This kernel has the following 
supermatrix elements 

(^A)pg,rs = J dXidX2dx[dX2(f>p(x[)f^{Xi) 



X/)(Xl, X2; x'l , X2)0;(X2)0.5(X2) 

A [{qr\wee\ps} - ^{qr\WeAsp}] , 



(B6) 



where {qr\wee\ps} are the two-electron integrals. The su- 
permatrix representation of the interacting polarization prop- 
agator xa is then found from the Bethe-Salpeter equation 
[Eq. (A17)] written in the spin-orbital basis 



B- A* 



1 



(B7) 



where and B,i are the so-called orbital rotation Hessians 

+A [{ib\Wee\aj) - ^{ib\Wee\ja)] , (B8a) 



(B-i)ifl,;ij = ^ [{ab\wee\ij} - ^{ab\wce\ji}] ■ 



(B8b) 



We need to consider the Unear response non-Hermitian eigen- 
value equation 



B^ 

B* A* 



1 

-1 



(B9) 



whose solutions come in pairs: positive excitation ener- 
gies a)„j with eigenvectors (X„,^,Y„,^), and opposite (de- 
)excitation energies -cl>„j with eigenvectors (y* ,,X* ,). 
Choosing the normalization of the eigenvectors so that 
X^_^X,„,i - yI^^Y,„j = 6nm, the supermatrix in,)(w) can be 
expressed as the following spectral representation (where the 
sum is over eigenvectors with positive excitation energies) 



1 



CO — (jJn,A + '0+ \ Y„,^ 
1 / Y 



[Ka (BIO) 



(jj + w„_,) - /0+ \ 



'■n.A 



(BID 
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The fluctuation-dissipation theorem [Eq. (A24)] leads to the 
supermatrix representation of the correlation part of the two- 
particle density matrix P^j (using contour integration in the 
upper half of the complex plane) 



P, 



c,A 



X! ,Y*"'", X! ,X*"'" 



(B12) 



the simple contribution coming from DloCw) resulting from 
its diagonal form [Eqs. (B3)], and the correlation energy 
[Eq. (A6)] has the following expression in spin-orbital basis 

1 r' 

Ec = - \ dAy^{ps\w\qr){Pc,A)pq,rs 

= ^ r fl'i2^(<//7|H',.|fl;)(Y„,,)L(YM);i> 

-'O ia,jb n y 
+{ij\Wee\ab){\„^,l)*^{Xn,A)ib + (fl^lWef IO')(X„,^)*,(Y„„i);fc 

+{ai\wee\ib) [(X„,^)*,(X„,,i)j/7 - 6ij6ah\ |, (B13) 

where out of the integrals {ps\w\qr) associated with the 
general perturbation operator only the integrals of the type 
{ib\wee\aj) associated with the two-electron contribution of 
the perturbation operator survive because of the occupied- 
virtual/occupied-virtual structure of the two-particle density 
matrix. Using now real spin orbitals, the correlation energy 
can be simplified to 

= - dAy{ib\Wee\aj){V,^,)ia.jh, (B14) 

2 -JO t^Jb 

where 

{^c,A)ia,jb - ^ (X„,i + ^n,A)ia (Xn,i + Y„_^) - Sjjdah, 
n 

(B15) 

or, in matrix form, 

P,,, = 2 (X„„, + Y„„,) (X,,,, + Y„,,f - 1. (B 16) 

n 

Using the well-known fact that, if + and A^ - are 
positive definite, the non-Hermitian eigenvalue equation (B9) 
with real spin orbitals can be transformed into the following 
half-size symmetric eigenvalue equation 



(B17) 



where M,, = (A^ - B^)'^^ (A^ + B^) (A^ - B^)'^^ and with 
eigenvectors/,,,^ = -y/cZi;;7(A^ - B^)"'^^ (X,,,,) + Y,,,^), and 

using the spectral decomposition M^'^^ - 2„ Wj^J,Z„,^Zj|^^, 
the correlation two-particle density matrix P^j can be ex- 
pressed as 



P,,, = (A^ - B,y'^ M:"^ (A, - B,)'/2 _ 1. 



2. Expressions for spin-restricted closed-sliell calculations 

For spin-restricted closed-shell calculations, the eigenvec- 
tors {X„j,Yn,A) can be transformed into spin-singlet excita- 
tion/diexcitation vectors 

('x„.,),„ = ^ [(X„,,),.^„T + (X„,,),^„J , (B19a) 



CynAa = ^ [(Y„,,),,,, + (YnAlai] ^ (B19b) 

and spin-triplet excitation/diexcitation vectors 

('•"x„,,)„ = -1^ [(X„ - (X„,,),,„J , (B20a) 



e'°yn,A)ia - ^ [(Y,,,,),,,, - (Y„,,),,,J , (B20b) 



( ' X„,^),fl — 0^n,A)i-^al ' 
C' ^yn,A)ia - (Y„„l),>-r ' 
( ' X„,^),a = (Xn,/l),|a| > 
C'^yn.A)ia - (Y„,^);|„j_ , 



(B20c) 
(B20d) 
(B20e) 
(B20f) 



the indices /, a, j, b referring now to spatial orbitals. With this 
transformation, the linear response eigenvalue equation (B9) 
decouples into a singlet eigenvalue equation 



with the singlet orbital rotation Hessians 

+A [2{ib\wee\ai) - ^{ib\wee\ja)] , (B22a) 



('B,,).^_ = A [2{ab\w,e\ii) - ^{ab\wee\ji)] , (B22b) 
and three identical triplet eigenvalue equations 

with the triplet orbital rotation Hessians 

(^^A,)._ = {ea - ei)6ij6ab - A^(ib\wee\ja), (B24a) 



(B18) 



(B24b) 
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Performing the sums over spins in the correlation energy ex- 
pression of Eq. (B14), one gets, for real spatial orbitals, 

1 r' 

£, = - dAy{ib\wee\aj)('Pc,A)ia,jb, (B25) 

where remains only the contribution from the spin-singlet- 
adapted correlation two-particle density matrix ('Pc,,) - 
2<ri=T.i I.cT2=tl(Pc,A)ia,aa,,jo-2ba2, which Can be Calculated sim- 
ilarly as before 



where ^M, = ('A, - 'B,,)'^'(1A, + 'Ba){'A, - ^B,) 



'P.,^ = 2 



n 



(B26) 
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